load package and example data

library(lsea)
# load example lipidomic composition data generated by adding random noise to real data
data(comp_data)
# load group labels (A or B) for differential testing and enrichment analysis
data(labels)

CLR transformation

#transform composition data (samples x features) using centered log-ratio transformation
# note: this is not for normalized abundance data!!
clr_data <- lsea::clr.transform(comp_data)
# subsequent functions expect samples as columns
clr_mat <- data.matrix(t(clr_data))

Differential composition testing

# to perform a paired t-test (samples must be ordered by pair)
t_df <- lsea::two.group.row.test(clr_mat, labels, test = "t", paired = TRUE)
# to perform unpaired wilcox test
# note: wilcox test will throw warnings for every row with ties, so suppress those warnings
w_df <- suppressWarnings(lsea::two.group.row.test(clr_mat, labels, test = "w", paired = FALSE))
t-test results ordered by p-value
stat mean1 mean2 dm pvalue padj
TG 56:4-FA20:3 4.396256 1.4517529 -2.892099 4.343852 0.0001654 0.1476720
TG 49:3-FA18:3 3.827157 2.1155015 -2.453107 4.568608 0.0007326 0.2358636
TG 54:5-FA16:0 -3.664517 -1.9583616 1.941328 -3.899689 0.0011143 0.2358636
FA 20:3 3.635246 1.3916176 -1.973312 3.364929 0.0012011 0.2358636
TG 52:6-FA20:5 3.598200 2.5091510 -1.210132 3.719283 0.0013206 0.2358636
PC 18:2_20:2 -3.244180 -0.7073557 2.084001 -2.791356 0.0032283 0.4616506
Wilcoxon test results ordered by p-value
stat mean1 mean2 dm pvalue padj
PC 16:1_18:1 625 3.0007067 -0.2891511 3.289858 0.0000017 0.0015611
FA 16:0 148 -0.3826287 2.6278213 -3.010450 0.0001078 0.0358698
CE 18:1 153 -0.5924634 2.7080542 -3.300518 0.0001607 0.0358698
PE P-16:0/16:1 576 2.4939921 -0.2272526 2.721245 0.0001607 0.0358698
TG 52:6-FA20:5 551 2.5091510 -1.2101318 3.719283 0.0009932 0.1446534
TG 46:1-FA16:1 550 2.1788510 -0.6821682 2.861019 0.0010621 0.1446534

LSEA

# we can use the t-statistic to rank lipid species and perform enrichment analysis using the GSEA algorithm
# note: the returned Wilcoxon statistic is not directional
res <- lsea(t_df, rnk_name = "stat")
Top 5 positive results
pathway pval padj ES NES
1 TG_UFA_12-16 0.0050970 0.3403915 0.4918451 1.715305
2 UFA_12-16 0.0050970 0.3403915 0.4918451 1.715305
3 TG_UFA_16 0.0332024 0.9229341 0.5311323 1.558536
4 CE_SFA_22-26 0.0059928 0.3403915 0.9663300 1.487966
5 dhCer_MUFA_22-26 0.0557103 0.9229341 0.8349662 1.439585
Top 5 negative results
pathway pval padj ES NES
279 FA_12-16 0.0354196 0.9229341 -0.7234549 -1.558625
280 CE_PUFA_18 0.0056225 0.3403915 -0.9360713 -1.604818
281 HexCER_MUFA_22-26 0.0233447 0.9229341 -0.7464305 -1.608124
282 HexCER_22-26 0.0233447 0.9229341 -0.7464305 -1.608124
283 TG_SFA_18 0.0274334 0.9229341 -0.6814703 -1.613043
284 CE_17-20 0.0030432 0.3403915 -0.6900139 -1.849869

Lipid species structure annotation

# LSEA works by annotating the provided lipid species based on their LIPIDMAPS-style name
lipid_anno <- lsea::annotate.lipid.species(rownames(t_df))
Species Class Category Total.Carbons Longest.Tail Total.DBs Saturation Chain
TG 56:4-FA20:3 TG.56.4.FA20.3 TG Glycerolipid 56 20 4 PUFA LCFA
TG 49:3-FA18:3 TG.49.3.FA18.3 TG Glycerolipid 49 18 3 PUFA LCFA
TG 54:5-FA16:0 TG.54.5.FA16.0 TG Glycerolipid 54 16 5 PUFA LCFA
FA 20:3 FA.20.3 FA Fatty.Acyl 20 20 3 PUFA LCFA
TG 52:6-FA20:5 TG.52.6.FA20.5 TG Glycerolipid 52 20 6 PUFA LCFA
PC 18:2_20:2 PC.18.2.20.2 PC Glycerophospholipid 38 20 4 PUFA LCFA
TG 52:4-FA18:2 TG.52.4.FA18.2 TG Glycerolipid 52 18 4 PUFA LCFA
CE 17:0 CE.17.0 CE Sterol 17 17 0 SFA LCFA
PE P-16:0/16:1 PE.P.16.0.16.1 PE.P Ether 32 16 1 MUFA LCFA
TG 56:6-FA22:4 TG.56.6.FA22.4 TG Glycerolipid 56 22 6 PUFA VLCFA
CE 18:1 CE.18.1 CE Sterol 18 18 1 MUFA LCFA
CE 18:4 CE.18.4 CE Sterol 18 18 4 PUFA LCFA
TG 49:2-FA17:0 TG.49.2.FA17.0 TG Glycerolipid 49 17 2 UFA LCFA
FA 16:1 FA.16.1 FA Fatty.Acyl 16 16 1 MUFA LCFA
TG 58:9-FA18:1 TG.58.9.FA18.1 TG Glycerolipid 58 18 9 PUFA LCFA

Visualizing differential lipid species enrichment

Plot the differential lipid species by longest tail length, number of double bonds, and class
# for the purposes of this example we will raise the adjusted p-value threshold
# since there are not many significant differences with this simulated data
lsea::structure.enrichment.plot(
  de_tbl = t_df, anno_tbl = lipid_anno, group_names = c("A", "B"),
  p_thresh = 0.96, color_pal = c("dodgerblue", "firebrick"), size_range = c(0.25, 2.5),
)

# or specify a Lipid Class
lsea::structure.enrichment.plot(
  de_tbl = t_df, anno_tbl = lipid_anno, group_names = c("A", "B"),
  p_thresh = 0.96, color_pal = c("dodgerblue", "firebrick"), size_range = c(0.25, 2.5), class = "TG"
)